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Abstract 

We study the transport equation describing a dense system of gluons, in the 
small scattering angle approximation, taking into account medium-generated 
effective masses of the gluons. We focus on the case of overpopulated systems 
that are driven to Bose-Einstein condensation on their way to thermalization. 
The presence of a mass modifies the dispersion relation of the gluon, as compared 
to the massless case, but it is shown that this does not change qualitatively the 
scaling behavior in the vicinity of the onset. 


1. Introduction 

In previous papers Di it was argued that a dense system of gluons such 
as those created in the early stages of an ultra-relativistic heavy ion collision, 
could be driven to Bose-Einstein condensation, as the system evolves towards 
thermal equilibrium. This was inferred from a detailed study of the kinetic 
equation that takes into account 2 to 2 scattering, in the small scattering angle 
approximation. Overpopulation means that the dimensionless number n/e 1 * 3 / 4 , 
where n is the number density and e the energy density, exceeds its value in 
equilibrium. An overpopulated system has too many gluons, relative to its 
total energy, to be accommodated in a Bose-Einstein distribution, and thermal 
equilibrium requires the formation of a condensate. 

Of course, such a condensate will develop provided the approach to thermal 
equilibrium proceeds with conservation of both energy and particle number. 
While energy is certainly conserved, inelastic processes of various kinds may 
change the number of gluons (see e.g. Enafl eventually preventing the forma¬ 
tion of a condensate in the true equilibrium state. However, particle number 
may be approximately conserved during much of the evolution, and this could 


1 Note that quark production, although it decreases the number of gluons, does not neces¬ 

sarily hinds the formation of a condensate]!]. 
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be enough to approach condensation. Indeed, transport calculations indicate 
that the amplification of soft modes is a very rapid process, that a chemical 
potential is indeed dynamically generated and that the onset for condensation 
can be reached on short time scales. This is confirmed by calculations using the 
small angle approximation [2], as well as complete solution of the Boltzmann 
equation urn- There are also indications that inelastic processes could accel¬ 
erate the amplification of soft modes |)j, while the authors of Refs. [5] and jJO 
seemingly reach a different conclusion. 

Clearly the analysis of inelastic scatterings requires further work, but this is 
beyond the scope of this paper. Motivation for studying the possibility for gluons 
to condense comes of course from the desire to better understand how matter 
produced in a high energy nucleus-nucleus collision evolves towards local thermal 
equilibrium (see (TT)1 ITT] for a recent review). But, as we already emphasized in 
[2], the general issue of the dynamical formation of a condensate is an interesting 
problem in itself. It is of relevance in the context of cosmology (see e.g. HU), 
or cold atom physics (see e.g. m)- It has be studied using kinetic theory, 
or classical field simulations (see e.g. ns ns)- In the context of Quantum 
Chromodynamics, the nature of the condensate remains an interesting, but 
unsolved, question (for a recent study in a related context, see ll6l). 

Our goal in this paper is to pursue our general study of the phenomenon 
within kinetic theory, using a transport equation that incorporates properly 
the effects of Bose statistics. The fact that the interactions are long range 
interactions validate the use of the small angle approximation which reduces 
the transport equation to a Fokker-Planck equation, much easier to solve than 
the Boltzmann equation, thereby providing more direct analytical insight. This 
paper, as well as a companion paper, addresses issues that were not discussed 
in Ref. [2] which is limited to the study of the onset of condensation. We want 
to extend our work so as to be able to obtain a complete dynamical description 
of the approach to equilibrium including the formation of a condensate. In 
order to do so, we need to attribute finite masses to the gluons. Such masses 
are automatically generated by the coupling to thermal fluctuations, and the 
proper transport equations that incorporate such self-energy corrections could 
be derived from first principles. However, for the purpose of the present study, 
it is sufficient to just give the gluons a mass, and correct appropriately the 
scattering matrix element. In fact two masses will be introduced. The screening 
mass niD regulates the infrared behavior of the collision kernel. The other mass, 
m, modifies the dispersion relation, and one of the issues that we want to study 
is how this modification changes the onset of Bose-Einstein condensation. We 
shall see that, in fact, qualitatively it does not. Finally, the role of the mass m 
is to allow a clear definition of the equations that describe the evolution of the 
system beyond the onset, that is, in the presence of the condensate. This will 
be discussed in a companion paper m- 

The outline of this paper is as follows. In the next section, we re-derive 
the approximate form of the transport equation in the small scattering angle 
approximation, paying attention to the presence of finite masses, and empha¬ 
sizing the differences with the massless case discussed in [2j. In section 3, we 
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present results of numerical solutions of the transport equation, illustrating the 
role of the mass on generic features of thermalization in both the underpopu¬ 
lated and overpopulated situations. In section 4, we focus on the critical regime 
that accompanies the onset of condensation. We show in particular that the 
change of the dispersion relation, from ultra-relativistic in the massless case, 
to non-relativistic in the massive case, does not change qualitatively the scaling 
regime. Details on the analytic calculation, as well as on the numerical solution, 
are given in two appendices. 


2. Derivation of the Transport Equation in the Massive Case 

In this Section, we derive the transport equation, under the approximation 
of small scattering angle, taking into account medium-induced effective masses 
for both the colliding gluons and the exchange gluons. 


2.1. Generalities 

As in our previous papers IDE], we assume a transport equation for the 
single particle distribution of the following form 

T>tfi = C[f} 

1 /' d 3 p 2 d 3 p 3 d 3 p 4 1 2 

2 J (2tt) 3 2E-2 (2tt) 3 2E 3 (2tt) 3 2E 4 2E ± 1 12 ^ 34 1 

x (2tt) 4 6(P 1 +P 2 -P 3 - Pi) 

X {/ 3 /4(l + /l)(l + /2Wl/2(l + / 3 )(l + /4)} (1) 


where 

V t = d t + v i-V, (2) 

and the factor 1/2 in front of the integral is a symmetry factor. Summation over 
color and polarization is performed on the gluons 2,3,4; average over color and 
polarization is performed for gluon 1. The distribution function / is a scalar 
object (i.e., independent of color and spin): 


f(x,p) 


(2tt) 3 dN 
2 (A/? — 1) d 3 xd 3 p' 


( 3 ) 


where N denotes the total number of gluons. In other words, / denotes the num¬ 
ber of gluons of a given spin and color in the phase-space element d 3 xd 3 p/ (27t) 3 . 
We consider a uniform system, so that / is in fact independent of x. Also, in 
this paper, we consider a non-expanding system, so that T>t = 9*. Finally, Pi 
denotes the four-momentum of particle i, Pi = (£’ p .,p i ). 

In this paper, contrary to our previous work, we assume that the gluons 
carry a small mass, arising from their interactions with the medium. This is 
a crude approximation to the true self-energy, but our main goal here is not 
to give a quantitative description of the phenomenon, but rather to study how 
the onset of Bose-Einstein condensation is affected by such a mass. We believe 
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that a more complete, but much more difficult, treatment would not change our 
main conclusions. The major change brought by the presence of the mass is 
that of the dispersion relation of the modes. Thus the energy of a particle with 
momentum p is E p = \J p 2 + m 2 . 

It is convenient to express the final momenta in terms of the initial ones and 
of the momentum q transferred in the collision. We set 

P 3 = Pi + 9, Pi = P 2 ~ 9 ( 4 ) 

One can then perform the integrations over p 3 and p 4 using the three-momentum 
delta function. We get 


C[f]= j f w(p 1 ,p 2l q){f 12 ->34}, 

Jp 2 J q 


with 


w{p 1 ,P 2 ,q) 


nS(E 1 + E 2 - E 3 - E a ) 
1 6Ei E 2 E% E 4 . 


I AI 12 —>34 Pi 


and 


{/l2^ 34 } = / 3 / 4 ( 1 + /l)(l + f 2 ) ~ /l/ 2 (l + / 3 )(1 + / 4 ), 


(5) 

( 6 ) 

(7) 


where /) is a shorthand notation for f p . The quantity w(p 1 ,p 2 , q) may be inter¬ 
preted as the rate of collisions of particle 1 with particle 2, in which momentum 
q is transferred to particle 1, its momentum p l becoming p 3 + q. Note that the 
symmetry factor is included in the definition of w. The (dimensionless) matrix 
element AIi 2 ->. 34 will be discussed shortly. It is understood in the expressions 
above that the momenta p 3 and p 4 are expressed in terms of p x , p 2 , q according 
to Q. Also, we used the shorthand for momentum integration 


/ 


d 3 p 

(2tt) 3 ' 


( 8 ) 


Note that the symmetries of the matrix element (see below) entail the property 


™(Pi,P 2 ,-<?) = w(p 1 ,p 2 ,q) = w(p 2 ,p 1: -q). 


(9) 


2.2. The small scattering angle approximation and the Fokker-Planck equation 
Under the small scattering angle approximation, the momenta of incident 
particles get changed very little during each collision, and the kinetic equation 
can be approximated by a Fokker-Planck equation in momentum space [15]. 
Following a standard procedure m, we write the collision integral as 

C[f\ = ~V ■J = - d ^. (10) 

vPi 

Note the gradient V in the above, and for the rest of this paper, is the momen¬ 
tum space gradient i.e. V = V p . In order to estimate , the component of the 
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current (of particles 1) in the direction i (with i = 1,2,3), we count the number 
of particles that, as a result of collisions during the interval d t, cross a surface 
element orthogonal to the direction i and located at a particular value p of p 1 . 
An elementary analysis yields 

Ji = f [ [ dp i w(p 1 ,p 2 ,q) 

JQ, Qi>0 J p 2 Jpi—qi 

x |/pi fp2 (-*■ + /pi+q) (1 + /p 2 -q) — (1 + /pi) (1 + /p 2 ) fpi+q /p2-q|' 5 

(11) 

where the q- integration is restricted to positive components qi. In the small 
angle approximation the combination of statistical factors simplifies into 

,fp 1 fp 2 (1 + fpi+q ) (1 + /p 2 ”<j) — (1 + fpi ) (1 + /p 2 ) fpi+q fp 2 -q 
*q-[h Pl (Vf)p 2 -h P2 (Vf) Pl ] +0(q 2 ), (12) 

where we have introduced the notation h p = / p (l + f p ). In addition, we notice 
that, in this approximation, the energy conservation implies 

0 = q ■ - q ■ v 2 + 0{q 2 ). (13) 

By taking all these together we obtain the following leading order expression 
for the momentum flux: 

3 l = \\ I Q l w{p 1 ,p 2 ,q)q • [V( v /)p 2 -V^iOpJ • ( M ) 

Z dq Jp 2 

Note that in the above equation we have relaxed the constraint q.j > 0 on the 
q-integration, dividing the result by 2 (using the fact that the integrand is even 
in q, see Eqs. <§)• We can rewrite the current as follows 

T=\ B^( Pl ,p 2 ) [h Pl (V\f) P2 - hp 2 (V j f) Pl ] , (15) 

with the (dimensionless) angular tensor 

= \ I Q l Q J w(p 1 ,p 2 ,q). (16) 


2.3. Conservation laws 

The particle number conservation is obvious due to the structure of the 
collision term as the divergence of a current: 

f C[f p J=0 (17) 

J Pl 
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To prove the energy conservation requires a little more effort 


[ E Pl C[f Pl ] — f E Pi V Pl ■ J(Pi) 

J Pl J Pl 

= - / V Pl • [E Pl J( Pl )] + f (V Pl E Pi ) • J{ Pl ) 

JP 1 -'Pi 

= / 

= o / / / w (Pi>P2;g) 


Pi J P 2 J 9 

x(wi ■ q) q ■ [V( V /)p2 - ^(V/)pi] 


o w{ P i, P 2\q) 

Z -'Pi ^P2 -'Q 

x [(wi ■ 9) h Pl (q • V/) P2 - (u 2 • q) /i Pa (q ■ V/) Pl ] , 

0 (18) 


where in the last steps we have used Eq.(13) and Eqs. <©• 


The matrix element 

The matrix element for (in vacuum) gluon-gluon scattering (1 + 2 
reads (spin and color averaged for 1, and summed for 2,3,4) 


M |"= 128t T z aiN z c 


tu su t s 
° ~2 72 77 


t 2 u* 

with s,t,u the standard Mandelstam variables: 
s = (-Pi + -P2) 2 , t = (-Pi — P3) 2 ) 


1287r 2 a 2iv2 = 72 g A 


= (Pi - P 4 


3 + 4) 

(19) 

( 20 ) 


In the small scattering angle approximation, the dominant contributions 
come from the kinematic regions t « 0 and u ~ 0. With the vacuum matrix 
element and in massless case, one would have the following approximation: 


I M | 2 = 72/ 


t u 

S U 

t s 


su t s 

3 2 

S 

t 2 

u 2 

-72 q 4 

t 2 u 2 



( 21 ) 


where in the last step we have used the fact that the two contributions are equal, 
and, for massless particles, u = — (s + 1) m — s (for t ~ 0). 

Now we consider the modifications of the matrix elements that need to be 
taken into account when the gluons are massive. There are two distinct physical 
effects. One is the screening of the t-channel singularity. This should be taken 
care of by separating the transverse and the longitudinal channels, and including 
the proper polarization tensors in the exchange gluons. In this paper we simply 
modify the denominator by substituting t -+ t — m 2 D , with mu a screening mass 
whose main role here is that of a regulator. The other physical effect is coming 
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from self-energy corrections on the external lines. As already mentioned, we 
simply take these into account by giving the gluon a small thermal mass m, 
allowing m to differ from mr>. In summary, we replace the matrix element 
derived for massless particles by 

| M | 2 —>■ 144g 4 g2 0 « 2E 1 E 2 (1 - Vl ■ v 2 ), (22) 

[t m D ) 

where v r = is the velocity of particle i, and since we are interested in the 

small m limit, we dropped a term ~ m 2 in the expression of s. With that, we 
then obtain, in the leading order of the small scattering angle approximation 

w = 36ng 5{q-v 1 -q-v 2 )-^ - 5 -(23 

(ur — — m D) 

where w = v± ■ q. 

With this matrix element, the angular tensor takes the form 

Bij (Pi,P2 ) = 18 V / S(vi -q~v 2 • q)q'q> ”2 )2 ■ ( 24 ) 


2.5. The isotropic case 

We assume that the distribution function is a function of the energy, f p = 
f(E p ) with E p = \J p 2 + m 2 . In such case we have the following relation 

Vf p = vf'(Ep), /'(£) = (V 

with the velocity given by v = p/E p = V p E p . Then we have 

h Pl (V J f ) P2 - h P 2 (V j f) Pl = v{ h l/' - v{h 2 f[, (26) 


where we introduced the simplified notation _/) = f(Ei) = f p . (and similarly for 
hi) that will be used throughout the paper. 

In this isotropic case, the calculation of the angular tensor (24) simplifies. 
In particular, the current J(p) is a vector aligned with the direction of p, that 
is (with pi = | p ± \ and p x = p 1 /pi) 


J(Pi)=PiJ(pi)- (27) 

It follows that the kinetic equation can be written as 

Vtfi = ~\d Pl {pU(pi)} . (28) 

Pi 

The calculation presented in Appendix A yields 


J (pi) = 36na 2 s f {h x f 2 - h 2 f{)Z(v 1 , v 2 ), (29) 

■>P2 

where the explicit expressions of the dimensionless function Z(v\,v 2 ) is given 
in Appendix A. 
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2.6. The massless limit 

It is shown in |Appendix A~T| that the current at small momentum has the 
following structure 


— ► 0 ) ~ 2>d>TTa 2 C 


I a (pi)d Pl fi + h(pi)hi 

£/i 


where £ is a positive constant, and (see Eq. ( A.25[ >) 

f'l 

—, ibW 1) = / 

'p 2 


IM = hM = f 

m J P2 v 2 m J P2 v 2 m J p 


\d p f 2 . 

< P2 


(30) 


(31) 


This structure is identical to that obtained in the massless limit [2], with the 
integrals I a and If, given by 


la = [ /(p)(l + m, h= f 

J r> Jr) 


2/(p) 


(32) 


We shall end this section by making a general comment on the Fokker-Planck 
equation in this massless limit, in order to clarify the interpretation of the two 
competing terms that are present in the current. For simplicity we focus on the 
isotropic case so that the kinetic equation is of the forir0 


= i^2 d P { P 2 l^dpfip) + hf{ 1 + /)]}■ 
Note that energy conservation, d t J p pf(p) = 0, entails 


df 

dp 


= h 


[ /(! + /), 

Jp 


(33) 


(34) 


a relation that is obviously satisfied given the definitions (32) of I a and lb- 
To gain insight into the physical meaning of the two terms on the r.h.s. of 
this equation, let us multiply both sides by p 2 and integrate over momentum. 
We obtain 

d t [ P 2 fip)= [ P 2 \d p {p 2 [I a d p f(p)+I b f(l + /)]}. (35) 

J p J n V 


Using integration by part for the r.h.s., we eventually obtain 
d t {p 2 ) = 6 I a (n) - 2 I b {p), 

where 


(36) 


(p)= P f(p), (n) = / f(p), (p) = / pf(l + f). 


(37) 


2 We absorb here the constant factor 367 tq 2 £ into the redefinition of time (see Eq. < |38| ) 
below). 










The equation (36) may be interpreted as follows. The first term on the r.h.s., 
6 I a (n), represents the diffusion of particles in momentum space that results 
from multiple small-angle scatterings. If this first term was the only one, the 
system would diffuse indefinitely, with ( p 2 ) ~ I a (n)t at late time. The second 
term on the r.h.s., —2 h(p), represents a drag force, with the characteristic of 
a friction. It is negative as long as there is any nonzero momentum in the 
system and therefore opposes the diffusive contribution, which is generically 
positive. In the absence of diffusion, this term would cause the magnitude of 
the momentum to decrease continuously, eventually bringing all particles to a 
state of zero momentum. This drag term is the agent that drives the system 
towards condensation. Finally, with both terms present on the r.h.s of the 
equation, thermal equilibrium can be reached, and the system approaches the 
fixed point corresponding to a Bose-Einstein distribution, with T = I a /Ib- 


3. Numerical solutions to the transport equation 


We turn now to the discussion of results obtained by solving numerically the 
Fokker-Plank equation for relevant cases. Details on the numerical procedure 
that we used can be found in |Appcndix B| We shall write the transport equation 
in terms of dimensionless quantities. In order to do so, we express all momenta 
(and masses, energies, temperature) in units of Q s , and the time in units of 
l/Q s . In fact, we also absorb into the dimensionless time a numerical factor, 
setting 


t 



■ 

18of 


(38) 


This is the factor that should be used if one wants to relate the time r of the 
simulation to the physical time t. [^] 

In the dimensionless variables, the transport equation reads 


Z>r/(P) = ~p2 d P {P 2 * J(.P)} > 


J(Pl) 



h2f'i)Z(v l ,v 2 ). 


(39) 


Note that since from now on we shall be dealing mostly with dimensionless vari¬ 
ables, we keep the same names for the dimensionful and dimensionless momenta. 


The initial condition is chosen be of the Glasma type DP, i- e - f ( E ) = 
/o(l + e 10 ^" 1 - 5 ))- 1 , with E = \fp^ + to 2 . For this family of initial condi¬ 
tions, there is a critical value of /o, that we call /q, above which the system 
becomes overpopulated. When /o = /g, the equilibrium distribution is a Bose- 
Einstein distribution with a maximal chemical potential /r = m. The value of 


3 Note that this factor d iffer s by a factor 2tt 2 C from the “natural” factor that appears for 

instance in the expression l29k of the current. Note that tt/(18o 2 ) r, 0.7 so that t r; t/Q s . 
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Figure 1: (Color online.) Left: Critical /o as a function of the gluon mass m. Points above 
the curve correspond to overpopulated initial conditions, below the curve to underpopulated 
initial conditions. Right: The evolution of /(p ~ 0) for different initial occupations /o (with 
m = 0.3): the lower (black) curve is for underpopulated initial condition, and the other curves 
are for overpopulated initial conditions. 


/o depends on the mass m of the gluons, and its variation with m is illustrated 
in Fig. [l] It is seen that f§ increases with m: the value of /o required to reach 
the onset of condensation is larger for massive particles than for massless ones. 
The right hand panel of Fig. [l] illustrates the behavior of f(p = 0) in the two 
generic situations of underpopulation where /(0) reaches a constant value and 
that of overpopulation where /(0) diverges at a finite time r = r c . 

Note that in all cases, the infrared part of the distribution function evolves 
rapidly towards an approximate classical thermal distribution function, 



(40) 


with an effective temperature T* and effective chemical potential /i* that can 
be determined numerically. The evolution of T* and n* with time depends 
on whether the system is over or underpopulated, as we shall discuss in the 
following. 

3.1. Underpopulated case 

Starting from the initial distribution corresponding to underpopulation, we 
verify that the system evolves to the expected equilibrium state. The approach 
to thermalization is illustrated in Figj2] which displays the evolution of the dis¬ 
tribution function f(p) and that of the current J[p). The current changes its 
sign around E ~ 1.3, the point which separates the effects of the two compet¬ 
ing components of the current. For larger values of E, the current is diffusive 
and positive. For smaller values of E, the current is dominated by its drag 
component, and is negative. As time progresses these two components of the 
current move particles in momentum space, the drag term pushing particles to¬ 
wards small momenta, the diffusion term smoothening the distribution at large 
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Figure 2: Left: Distribution function f(E) at various times r =0.05, 0.8, 1.4, 6.25 (from 
bottom to top at E = 0.4). Right: Current J at r =0.05, 0.8, 1.4, 6.25 (from bottom to top 
at E — 0.4). Underpopulated initial condition, with m = 0.3 and /o = 0.1. 


momenta. The resulting distribution gradually evolves towards the equilibrium 
distribution, as indicated in the left panel of Fig. [2j while the current diminishes, 
and eventually vanishes (which takes place approximately for the largest time 
considered in the plot). The evolution with time of the effective parameters that 
characterize the infrared part of the distribution function are displayed in Fig. [3] 
This figure clearly demonstrates that the system tliermalizes as expected. 




1 2 3 4 5 6 

X 


Figure 3: (Color online.) The evolution with time of the effective local temperature T* 
(left) and chemical potential fi* (right) for an underpopulated initial condition (m = 0.3 and 
fo = 0.1). At late time, both quantities approach their expected equilibrium values, indicated 
by the horizontal (red) dotted lines. 


3.2. Overpopulated cases 

We turn now to the situation with overpopulation. We expect, from our 
previous work, that the system will approach the onset of condensation in a 
finite time r = r c . We have studied the evolution with different masses(m = 
0.1, 0.3, 0.5, 0.7). In the numerical solution, the time step is self-adaptive i.e. 
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becoming smaller when the evolution becomes faster, which allows us to evolve 
the system in each case very close to the onset. The evolution is stopped when 
/ > 800 for the smallest momentum grid point, which roughly corresponds to 
the situation when the difference between the local chemical potential /z* and 
the mass becomes less than the mesh size. That is, the evolution is stopped 
very near the onset of condensation where /z = to. An illustration of the energy 
dependence of the distribution near the onset is displayed in Fig. [4] 



Figure 4: The inverse of the distribution as a function of energy, 1/f versus E , for m = 0.3, 
and r = 0.101 corresponding to the onset of condensation. 


The Table [l] gives some typical values of the parameters /g and r c for differ¬ 
ent initial conditions, and different values of the mass m and the Debye mass 
mu. One notices that, for fixed masses (e.g. m=0.3) the onset time decreases 
with increasing / 0 , which is as expected: the larger the overpopulation, the 
shorter time it takes to reach the onset of BEC. For fixed to and /o an increase 
of the Debye mass leads to an increase of r c (compare for instance the lines 
corresponding to m = 0.7 and /o = 1). This is because an increase of tod 
reduces the scattering rate, which slows down the evolution. A larger increase 
of t c results from the increase of m at fixed /o and mjj. The latter phenomenon 
is in line with what we observed earlier, namely that /g increases with to. 

Different from the underpopulated case, the evolution of the flow of particles 
is quite different in the overpopulated case, as shown in Fig(5] . The negative 
infrared flux keeps increasing in a roughly self-similar manner. The local chem¬ 
ical potential /z* approaches /z c = m at the onset. We have extracted the /z* 
as a function of time for a variety of mass values as well as initial occupation 
values: see Fig(6] In all cases, we have found that the evolution of /z* close to 
the onset point can be fitted with power law /z = to — A (r c — t)A Furthermore 
in all cases, we’ve found the optimal exponent is about one, i.e. i] — 1. This 
dynamical behavior in the present massive case is essentially the same as that 
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m 

m D 

/o 

fo 

T c 

0.1 

0.1 

1.0 

0.271 

0.029 

0.3 

0.3 

1.0 

0.508 

0.101 

0.5 

0.5 

1.0 

0.747 

0.259 

0.7 

0.7 

1.0 

0.969 

0.663 

0.3 

0.3 

1.2 

0.508 

0.065 

0.3 

0.3 

1.4 

0.508 

0.046 

0.3 

0.3 

1.6 

0.508 

0.034 

0.3 

0.3 

1.8 

0.508 

0.026 

0.3 

0.3 

2.0 

0.508 

0.021 

0.1 

0.7 

1.0 

0.271 

0.095 

0.7 

0.1 

1.0 

0.969 

0.242 


Table 1: Critical occupation factor f§ and onset time scale r c for different combinations of 
parameters. 




Figure 5: Color online. Distribution function / (left) and current (right), the curves 
corresponding, from bottom to top, to increasing values of r : r =0.04(black), 0.0674(red), 
0.0871 (green), 0.0927(blue), 0.0961 (yellow), 0.0977(pink) with m = 0.3 and / 0 = 1.0 (over- 
populated initial condition). 


found for the massless case in 5]. We give a more complete discussion of this 
regime in the next section. 

4. Critical scaling analysis 

In this section we examine how the scaling behavior that holds near the 
onset of BEC is affected by the gluon mass. A major modification, with respect 
to the massless case, concerns the dispersion relation near the onset. When 
m 0, condensation occurs when m = ji, and the dispersion relation becomes 
non relativistic, E p — p, = \/p 2 + to 2 —to ss p 2 / (2m), which differs from the ultra 
relativistic dispersion relation of the massless case, E p = p , where condensation 
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Figure 6: Evolution of IR local chemical potential //' toward onset /i* —> m for different 
choices of /o =1.0, 1.2, 1.4, 1.6, 1.8(from left to right) with the same mass m = 0.3. Right: 
Evolution of IR local chemical potential y* toward onset y* — > m for different choices of 
external mass m =0.7, 0.5, 0.3, 0.1(from left to right) with the same fa = 1. 


occurs at /it = 0. This change in the dispersion relation modifies the singularity 
near the onset, but, as we shall see, this does not affect in a major way the 
onset critical behavior. The foregoing analysis follows closely that presented in 
Ref. PJ. 


f.l. Scaling behavior of the current at small momentum 

As we have seen in solving numerically the Fokker-Planck equation, for 
generic initial conditions, the distribution function at small momentum evolves 
rapidly towards the approximate equilibrium distribution: 


f(j>) 


E{p ) - p* 


2 mT* _ 2 mT* 
p 2 + 2 mSp p 2 + A 2 ’ 


(41) 


where we have set Sp = m — p* and A = y/2m6p. As the system approaches 
condensation, A —> 0. To analyze the detailed behavior of the properties of the 
system when A —> 0, it is convenient to calculate the time dependence of the 
number of particles inside a small sphere of radius po centered at the origin. 


In doing this calculation, we assume that f(p) keeps the form (41), with time 
dependent parameters T* and Sp. In particular, /(0) remains finite as long as 
Sp, ^ 0. A simple calculation then yields 

rpo 

dr dpp 2 f(p) = 2mA(d T T*)h 1 (y)+2mT*(d T A)h 2 (y), (42) 

J o 


where y = Po/A , and the two scaling functions are 


hi(y) = y — ArcTan(j/) (43) 

h 2{y) = * f 2 - ArcTan(y). (44) 

1 + y 2 
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By using Eq. (281, one can relate the left hand side of Eq. (42) to the current 
at po , and obtain 


- J(po) = 


dr fo° dpp 2 f(p) 2m(d T T*) hi(y) , 2mT*{d T A) h 2 (y) 


P 2 o 


y- 


A 2 


r 


(45) 


This equation provides interesting constraints on the small p behavior of J(jp). 
In the limit y —► 0 or p 0 <C A, one gets 


_, . 2 m _ 

~i7(Po) — O t 


rjn* 

A 2 


= Po x -d T 


Ji* 

5/i 


= y3 r /(0). (46) 


The current is linear in p^ 1 with a slope proportional to d T f( 0), very much like 
in the massless case [3]. On the other hand, in the limit!/ > 1 or A <Cp 0 < T*, 
one obtains the following leading order result^] 

- J(Po) -X ~2 x (7rm)9 r (-T*A). (47) 

Po 


In this region, the current exhibits a singular behavior in 1/pg. The value 
(« A) of the momentum where the change of regime occurs decreases with 
time, while the absolute value of the current at the minimum becomes larger 
and larger, and eventually diverges at the onset. This behavior is illustrated in 
Fig.0 As was the case in the massless case, the current in the scaling regime is 
dominated by the second term in Eq. (451. One also finds that the dependence 
of to — p* on t c — t is linear. A discussion of the latter point is presented in 
Appendix |Appcndix A.4| 


5. Summary 

In this paper, we have extended the analysis carried out in Ref. [2] to the case 
where the gluons are given a small mass m. The main difference with respect 
to the massless case studied in |2j is the change of the dispersion relation of the 
gluon modes, from ultra-relativistic in the massless case to non-relativistic in the 
massive case. While this change in the dispersion relation leads to an enhanced 
singularity of the distribution function near the onset of Bose condensation, from 
1/p to 1/p 2 , we find that the critical regime that accompanies the approach to 
the onset is qualitatively unchanged. The Debye mass mn controls the infrared 
behavior of the collision kernel. It is kept constant, and independent of m, 
although in a more complete theory, both masses are related, and would depend 
on the temperature. However this is expected to have little impact at the onset, 
since the onset process is dominated by a critical regime where the actual value 
of the mass is irrelevant. We have checked in articular that letting the Debye 


4 The subleading contribution is ~ 1/p and comes from the first term in Eq. ( |45| ). It 
becomes the dominant contribution after onset, when Sp, = 0 and therefore (|47|) vanishes. 
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Figure 7: The current (in arbitrary units) as a function of momentum, at different time 
moments (earlier to later time from top to bottom) close to onset. Left panel: calculated 
according to the formula ( |45[ ), assuming a linear relation between //* — m and r c — r, and 
neglecting the time variation of T*; right panel: obtained from the numerical solution. The 
linear behavior at small p followed, as p increases, by a singular behavior in l/p a , with a > 0 
is clearly visible on both figures. 


mass adjust with temperature as one approaches BEC does not lead to any 
significant changes in the behavior of the system. Finally, the presence of the 
mass m allows us to define properly the equations that governs the dynamics 
beyond the onset, that is, in the presence of the condensate. This is discussed 
in a companion paper [El- 
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Appendix A. Derivation of the current J ( p) 


The component i of the current J7(p x ) in Eq. (15) can be written 

'■ d 3 q f d 3 p 2 


J\Pi) = 187 rg 


(2t tY J (271-)' 


M - h 2 fi) 


(A.l) 


x q l (q-v 2 ) 


(1 - ui • v 2 y 


[(«i • qf ~ q 2 - ™ 2 d\ 


■S(q-v 1 -q- v 2 ). 


16 

















Appendix A.l. Details of the angular integration 

In order to perform the angular integrals, we choose the following coordinate 
frames: We define the orientation of q in a frame where p 1 is along the z axis, 
and denote the corresponding angles by 9 and <fr. The orientation of p 2 , given 
p l and q , is defined by the angles 62 , <p 2 in a frame with q along z 2 and q,p 1 
spanning the x 2 -z 2 plane. We have therefore q v 1 = qiq cos 9, qv 2 = qv 2 cos 0 2 , 
and v 1 ■ v 2 = r'i , i’ 2 (cos#cos# 2 + sin#sin# 2 cos<^ 2 ). 

By integrating over <f> 2 one gets 

/ 0 27r d<j) 2(1 — Vi ■ v 2 ) 2 = d(j> 2 [1 — viv 2 (cos 9cos9 2 + sin0sin# 2 cos<^ 2 )] 2 

= 7r[2(l — v\v 2 cos 9 cos 9 2 ) 2 + w 2 u 2 (l — cos 2 0)(l — cos 2 6* 2 )] (A.2) 

Furthermore we rewrite the delta function as 


qd(q ■ Vi — q ■ v 2 ) = q6(q(v 1 cos 9 — v 2 cos 0 2 )) = S(v 1 cos 9 — v 2 cos9 2 ) (A.3) 
= 9(v 1 — v 2 )5[v\ cos 9 — v 2 cos9 2 ) + 9{v 2 — v\)5[v 2 cos 9 2 — v\ cos 9) 


and change variables, setting x\ = Vi cos 0, x 2 = v 2 cos 9 2 . One can then com¬ 
plete the angular integration. Note that, by symmetry, the current is aligned on 
the direction of p x , i.e., J[p\) = P\J(p \), with J{pi) a function of p\ = |p x | 
only, and p l = p 1 /p\.. We get 


J(pi) = J dp 2 p 2 2 (hif ' 2 - h 2 f[) 

with 

ViV 2 Z(vi,v 2 ,c q ) = - 


d<z Z(v 1 ,v 2 ,c q ) 


q Vi 


&X 2 X 2 


2(1 -x%) + {vi - x%)(v% - x%) 


= V 


j'V 1 

— Vl 
2 


~V2 [X\ - Cg) 

2 2(! - x\f + {vl - xl)(v% - x\) 


(A.4) 


9{v 1 - v 2 ) 


dx\Xi 


{x\ - CgY 


9{v 2 - vi) 


V — 6(1 — Cq) + (1 — Vl) + (1 — V 2 ) 


(j - Cqf (C q ~ vl)( Cq ~ vl) 

Cq - V 2 2 {c q - V 2 ) 

(! - Cqf _ {c q - VpjCg - Vp 
Cn 2Cn 


+ 6(1 — Cq) y/c^ArcTanhl 


\s/C Q 


-(l-u 2 )-(l -vl) 

(A.5) 


where c q = 1 + m 2 D /q 2 1 v = min(r)i, v 2 ), and we have used the following integral 
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formula to get the final results 


dx i 


3 C 


= Car — B x — 


B A 

Co - X 2 (c 0 - x 2 ) 2 
Ax 


(A.6) 




ArcTanh 


x 

V^oJ ’ 


2(x 2 - c 0 ) 

with A = 2(cq — l) 2 + (cq — f 2 )(co — v 2 ), B = 4 + v\ + v\ — 6cq, C = 1. 


Appendix A.2. Limiting cases 

There are a few limiting cases where the expression above simplifies. For 
instance, in the limit of vanishing screening mass mu —> 0, c q = 1 and 


V\V 2 Z = v 


+ (l-u 2 ) + (l-u 2 )- 


(l-ui 2 )(l-u 2 )i 


2(1 — v 2 ) 


(l-v 2 )(l-v 2 ) 


-(l-v 2 )-(l-v 2 ) 


ArcTanh(u). (A.7) 


In this limit, the remaining g-integration in the expression (A.4| of the current 
becomes simply the usual Coulomb logarithm, J dq/q —> C. 

In the limit where all effective masses vanish, i.e. m —> 0 and mu — > 0, c q 
and all the velocities become unity, V 1 V 2 Z —> 1 and we recover the expression 
of the massless case studied in Ref. [5J. 

In the limit where the colliding gluons become massless, i.e., m —>■ 0, but the 
Debye mass stays finite, all the velocities v±,V 2 ,v approach unity, and we get 


vi v 2 Z 




+v^i?(l 


Cq) 


3(1 ~ Cq) 
2 Cq 



ArcTanh 



(A.B) 


This expression is finite when q —» 0. To see that, recall that when q —>■ 0, 
c q ~ m 2 D /q 2 —» oo. As simple calculation then yields 


v x v 2 Z « Ac q 2 + 0(c q 3 ), 


(A.9) 


where 




(A.10) 


This result confirms that in the presence of a non-vanishing screening mass mu, 
the function Z is regular as q —> 0, and the g-integration in the expression (A.4| 
becomes infrared finite, as expected. 


18 





















Appendix A.3. Detailed derivation of the q-integration 

In the general case, in order to perforin the ^-integration in Eq. (A.4), we 
rewrite Viv 2 Z in Eq. (A.51 as follows: 

viv 2 Z(c q ,v 1 ,v 2 ) = 2 + E 2 (e q - v 2 ) + E 3 + y^ArcTanh ( ] 

c q V \ JC q J 


X [.£ 4(1 — Cq) + Eq — -b Eq\, 


(A.H) 


where 


E x = ^[2u' i (3v 2 -t;f-t^-4) + 2i;^+4], 

a = 1 


E, = 


15v 


E 3 = -(20v 2 - 3vl - 3v% - 12), 

£5 = ^[-2 v\vl - 4], E e = * [-2u 2 u 2 + 6u 2 + 6v% - 10]. (A.12) 

Then, we change integration variable from q to c q and obtain 

J ^-viV 2 Z(q,vi, v 2 ) = J 2 (i°- c \ v ^ Z i c q ,^ 1 ,^ 2 ) 

= £iGl(u,Cg) + E 1 G 1 (v,Cq) + E 2 G 2 (v,C q ) + E 3 G 3 (v,c q ) 
+£ 464 ( 1 ;, Cg) + EqGq(v, Cq) + EqGq(v, Cq). (A.13) 


The results of these integrals are listed below: 


2(l-Cg)(cg-i; 2 ) 


c a -v 


qtV'q 

2 1 


Gl(v,C q ) = J dc q 
G 2 {v,C q ) = J dc q 
G 3 (v,c,) = j ^ 

G 4 (v,c q ) = Jdcg^y/c^ ArcTanh f- 


ln(c 9 - v 2 ) - ln(c g - 1) 
2(v 2 - 1) 


2(1 - Cq) 


= ~2\ Cq + ( 1 ~ V ) ln (^<? x )] 


1 1, , 

■q on --T = ~o ln ( C 9 “ !) 


= - [2 (?J 2 ArcTanh ( - 


6 9 V+Cg 


1 


G 5 (v,c q ) = / dc q — ArcTanh —— = - 


' 2 c, 


1 


uln(c q — v 2 ) + 2 v +^ArcTanh I - 


+ V 3 In (c q - V 2 ) + C q v] 


(A. 14) 
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G 6 (v,c q ) = /^^l^ArcTanh 


= l ’{- : 

1 

4v 


-ArcTanh 


v 

./Cn 


In 


v 

,fCr, 


— - In ( 1 — 


,/Cq — 1 / V \ /I — vL/Ca\ ( V + V,/Cn\ 

21 ”<5#TT )AreT “ h hi) +Lb ( i+v ) +Li » (-- Gv 1 ) 


+Li 2 (^ 


V—V^Cq s , l + w/,/c ( 


-) + Li 2 (- 


1 — v 1+v 

+ ln(v - v/y/c^) ln(l - v/yfc, q ) - ln(u 2 - v 2 /c q ) ln(l - v) 

- ln(l - v 2 /c q ) ln(l + v) + ln(u + v/y/c^) ln(l + v/y/c^)] }. 

We define, with A an ultraviolet cutof0 
f A 

/ —viv 2 Z(c q ,vi,v 2 ) = L(ca,vi,v 2 ) - L(oo,Vi,v 2 ), 
Jo 9 


(A.15) 


(A.16) 


where the first argument of L denotes the value of c q (when q = A or q = 0, 
respectively). As we have already shown, there is no singularity at q —> 0 i.e. 
c q —^ oc when tod is finite, and we have the following finite result 


1 


L(oo,ui,u 2 ) = \ 20u 3 - I2v(v\v\ +2)] 


—3 A 


Li 2 


v — 1 


+ Li 2 


v + 1 


+ 2v — (2v + ln(l — v)) ln(u) 


(A.17) 


where A = —2v 2 v 2 + 6(vf + v 2 ) — 10. Note in the massless limit for external 
gluons, i.e. m —> 0, one has A —► 0 and the expression above greatly simplifies: 


L(c a ) - L(o o) = - - 


5c a + (3 - 5c A )v / CAArcCoth( v /cI) - - 


We write the final result for the current in the form 
18a 2 r 


J(P i) = 


dp 2 p\ Z(vi,v 2 ) [/ 11/2 - h 2 f[\ 


(A.18) 


(A.19) 


with 


Z{v i,u 2 ) = J 


d g Z(v 1 ,v 2 ,Cg) 
q vi 

2 n , \-l 


= {v(v 2 ) [L(ca,v-l,v 2 ) - L(oo,Vi,V 2 )} . 


(A.20) 


5 One expects A to be typically of the order of the temperature. The numerical calculations 
have been performed with a somewhat larger value, A = 4 Q s . 
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Appendix A A- The small momentum regime 

We now analyze the small momentum form of the current, J(j>\ —¥ 0). The 
leading order expression of Z in the limit v± —> 0 can be obtained after a lengthy 
but straightforward calculation. It reads 

v 3 

Z(v l ^t) 1 v 2 )^ (A.21) 

v(v 2 

where v = min(wi, U 2 ), and £ is a positive constant, 


£ = -3 


log(c A - 1) - log(c A ) H- 

CA 


CA = 1 + ^. (A.22) 


The small momentum current can therefore be written a^£] 

-J(pi) = 36tt a 2 s C [ [h 2 f[ - ft 1 / 2 ] ■ 

J P2 v x v 2 


(A.23) 


Since the numerical factor in front of the integral is to be absorbed in the 

(A.24) 


redefinition of the time scale (see Eq. (38)), we rewrite J{pi) as 

= Ia(Pi)fi+h(Pi)hi, 


where 


laiPl) = [ ft'2 ^ > Ib(Pl) = ~ [ 

J P2 l\V 2 J p 


/ /2-2- 
p 2 V l V * 


(A.25) 


These integrals reduce to the integrals I a and h of Eq. (32) in the massless 


limit (in this limit, the dependence on pi disappears and C becomes the usual 
Coulomb logarithm). 

We argue in the main text, and we have verified through numerical calcu¬ 
lations, that at small momentum and near the onset for BEC, the distribu¬ 
tion function f(p) can be well approximated by a Bose equilibrium distribution 
function We therefore write the distribution in the momentum range of interest 
(p«A = y/2m8p) as 


fip) = f* ( P ) + 8f(p), (A.26) 

where f* (p) is an equilibrium Bose distribution with temperature T* and chem¬ 
ical potential p *, and Sf(p) represents a small deviation from this equilibrium 
distribution. We assume that 8f(p) is a regular function of p at p = 0. 

In order to calculate the integrals I a b we need to pay attention to the fact 
that v = min(i>i, v 2 ), and divide the p 2 integration range appropriately, making 


6 It is not difficult to show that, with this approximate expression for the current, particle 
number as well as energy are conserved. 


21 












explicit the dependence on pi. We obtain 


27r 2 Iaipi) 


rp 1 „2 

dp 2P l^h 2 


dp2pl~h 2 

v 2 


2 I V 2 Vl 


dp 2 p 2 ~ ) h 2 + 

L ^2 


dp 2 p 2 2—h 2 

v 2 


np i / 2 \ /*00 

= / dp 2 p 2 [ ^ ) h’2 + — / dp2p 2 E 2 h2. (A.27) 

Vo VH P2/ m Jo 

In the first integral, we could set iq ss Pi/m and v 2 = p 2 l'rn, while in the 
second integral we need to keep the exact expression v 2 = p 2 /E 2 . Since, when 
Pi <C \/2m5p, h 2 in the first integral is nearly a constant, this first integral is 
of order p \, and is therefore negligible compared to the second one oc pi/m. We 
are then left with the integral given in Eq. (31). Proceeding in the same way 


for — J;,, one obtains an identical expression to that of I a with f 2 substituted to 
h 2 . We shall denote by I* and /£, the integrals I a and /{, calculated with the 
distribution /*, and write accordingly I a b = f* k + <57 a ,f>- We may then expand 


the current (A.24) as follows 


- J{pi) = SI a (f*y + 8I b h\ + r a {8fi)' + I{:8hi, 


(A.28) 


where Shi ~ Sfi(l + 2/j"), and we have used the fact that I* = T*I£ which 
follows immediately from the relation (/*)' = —(1 /T*)h*. By using again the 
same identity for the first two terms of the equation above, we can rewrite the 
current as 

J(Pi) = ~(6Ia - T*5I b )h\ - P a {5fi)' - 1^5h\ ■ (A.29) 

At this point, we note that the last two terms in the expression above can be 
neglected as pi —> 0. Indeed, 5f(jpi) is regular as pi —» 0, while h\ ~ (/i) 2 
diverges {Shi ~ Sfif* is subleading). A simple calculation yields 


Si a - SI b T* « ^ / dp 2 p 2 E 2 [Sf 2 { 1 + 2/*) + T*{5f 2 )']. 


(A.30) 


Thus, after dropping the last two terms in Eq. (A.291 one can rewrite J as 
follows 


= ^7( t )pi/ 2 (pi)i 


(A.31) 


where j{t) is a priori a regular function of time which can be expanded around 
the onset time r = r c , y(r) ~ y c + a{r c — r). By inserting this expression in 
Eq. (46) one gets the following equation for /(0) ~ f(pi) 


9rf Pl =1 (t), 


(A.32) 


whose solution reads 


/ p ! 1 ^) - - 7 c{t c - t) - -a{r c - r ) 2 


(A.33) 
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x X 


Figure A.8: / P1 1 (left) and 7 (right) as function of time with m =0.3. Red lines are linear 
fitting results. 


where we have used the expansion of 7 (r) and the condition / j 7 1 1 (t c ) ~ 0. We 
can check numerically that this is the correct behavior. In order to do so, we 
first fit the calculated function /^ t 1 (r) near r c with a linear form, 7 c (r c — r), 
and get the onset time t c ~ 0.1006 as well as 7 C ~ —0.6. Then we compare 
these values to those obtained from the fit of the quantity 3 (I a — hT*) / (p±T*) 
which exhibits a linear time dependence near r c , as expected. The two values 
agree perfectly, as can be seen in Fig |A. 8 | 

It is interesting to compare the present analysis to the analogous one pre¬ 
sented in Ref. [2|. In this case, pi/rn —> 1, and the combination I a — T*Ib 
vanishes. Both in the massive and the massless case, one finds that p* — m 
vanishes linearly with r c — r, but while this result follows from a simple argu¬ 
ment in the massive case, in the massless case this could only be determined 
numerically. 


Appendix B. Appendix B: Details of Numerics 

In this Appendix, we give some details on how we solve the Fokker-Planck 
equation. An efficient strategy is to solve its “integrated version”. Namely, 
instead of solving the differential equation for the distribution function directly, 
we evaluate the total number of particles in a thin shell dpp 2 f{p) ~ 

p 2 Apf{p) and examine its time evolution by integrating the transport equation 
over this moment window: 

p 2 Apd T f(p) = (p 2 J)p±Ap/2 = ^(p + Ap/2) - Tip - Ap/2) (B.l) 

Note on the right-hand side the kernel (taking the form of a full derivative 
~ V • J) will be integrated to give the difference of the flux T = p 2 J on the 
two surfaces of this shell at p ± Ap/2 respectively. 

Numerically we discretize the distribution on an equally spaced momentum 
grid pi = (i— 1/2)Ap, i = 1,..., 400, where A p = 0.01. So the flux is on the grid 
Pi = iAp, i = 0,1, ...,400. It is easy to see that T{p = 0) = 0. We further set 
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the flux to vanish at our momentum grid’s UV cutoff, i.e. SF{p = A = 4 Q s ) = 0, 
as the boundary condition that enforces exact particle number conservation. To 
solve this equation we use the implicit Gaussian scheme (which is a standard 
algorithm for this type of equation) as 

St 

Jt+St (Pi) ~ /r(Pi) p2/\p (^~Pi+Ap /2 [fr+fir} ~ ^~Pi — Ap/2 [/r+5r]) (B-2) 

where j c p [f T +sr] is the flux at momentum p evaluated with the distribution at 
t + St. The implicit scheme then involves numerically solving the above set 
of equations to extract the distribution at time r + St. The advantage of this 
implicit scheme, as is well know, is its robust numerical stability as compared 
with e.g. the explicit scheme of directly evolving the equation in time. During 
the whole evolution we have implemented an automated adjustment of the time 
step St to guarantee that the distribution function at the lowest grid point 
(which has the largest occupation f(p)) p = Ap/2 varies less than 5% at each 
time step forward. This allows rather accurate handling of the very infrared 
part of the evolution which is important for understanding the critical behavior 
in the onset of condensation. In the entire calculation the particle number 
conservation is exact while the energy conservation is maintained at the order 
of 10~ 3 variation or less. 


Reference 

References 

[1] J. -P. Blaizot, F. Gelis, J. -F. Liao, L. McLerran, and R. Venugopalan, 
Nucl. Phys. A 873, 68 (2012). arXiv:1107.5296 [hep-ph]]. 

[2] J. -P. Blaizot, J. Liao, and L. McLerran, Nucl. Phys. A 920, 58 (2013) 
arXiv:1305.2119 [hep-ph]]. 

[3] A. H. Mueller, A. I. Shoshi, and S. M. H. Wong, Phys. Lett. B632, 
257-260 (2006). hep-ph/0505164|; A. H. Mueller, A. I. Shoshi, and 
S. M. H. Wong, Eur. Phys. J. A29, 49-52 (2006). |hep-ph/0512045 . 
A. H. Mueller, A. I. Shoshi, and S. M. H. Wong, Nucl. Phys. B760, 145-165 
(2007). [hep-ph/0607136 . 

[4] J. P. Blaizot, B. Wu and L. Yan, Nucl. Phys. A 930, 139 (2014) 
arXiv: 1402.5049 [hep-ph]]. 

[5] F. Scardina, D. Perricone, S. Plumari, M. Ruggieri and V. Greco, Phys. 
Rev. C 90, no. 5, 054904 (2014) |arXiv:1408.1313 [nucl-th]]. 

[6] X. G. Huang and J. Liao, arXiv:1303.7214 [nucl-th]. 

[7] Z. Xu, K. Zhou, P. Zhuang and C. Greiner, arXiv:1410.5616 [hep-ph]. 


24 






[8] A. Kurkela and G. D. Moore, Phys. Rev. D 86, 056008 (2012) 
arXiv: 1207.1663 [hep-ph]]. 

[9] M. C. A. York, A. Kurkela, E. Lu and G. D. Moore, Phys. Rev. D 89, 
074036 (2014) |arXiv:1401.3751 [hep-ph]]. 

[10] J. Berges, J. P. Blaizot and F. Gelis, J. Phys. G 39, 085115 (2012) 
[arXiv: 1203.2042 [hep-ph]]. 

[11] X. G. Huang and J. Liao, Int. J. Mod. Phys. E 23, 1430003 (2014) 
arXiv:1402.5578 [nucl-th]]. 

[12] D. V. Semikoz and I. I. Tkachev, Phys. Rev. Lett. 74, 3093 (1995); Phys. 
Rev. D 55, 489 (1997). 

[13] C. W. Gardiner, M. D. Lee, R. J. Ballagh, M. J. Davis and P. Zoller, Phys. 
Rev. Lett. 81, 5266 (1998); C. W. Gardiner, P. Zoller, R. J. Ballagh and 
M. J. Davis, Phys. Rev. Lett. 79, 1793 (1997). 

[14] J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan, 
arXiv: 1303.5650 [hep-ph]. 

[15] T. Epelbaum and F. Gelis, Phys. Rev. Lett. Ill, 232301 (2013) 
arXiv:1307.2214 [hep-ph], arXiv:1307.2214 [hep-ph]]. 

[16] T. Gasenzer, L. McLerran, J. M. Pawlowski and D. Sexty, Nucl. Phys. A 
(2014) arXiv:1307.5301 [hep-ph]]. 

[17] J. -P. Blaizot, J. Liao, and L. McLerran, in preparation. 

[18] A. H. Mueller, Nucl. Phys. B 572, 227 (2000) |arXiv:hep-ph/9906322 ; 
Phys. Lett. B 475, 220 (2000) arXiv:hep-ph/9909388|. 

[19] E.M. Lifshitz and L.P. Pitaevskii, Physical Kinetics (Pergamon Press, New 
York, 1981). 


25 




